# data wrangling
library(tidyverse, warn.conflicts = FALSE, quietly = TRUE)
library(tidyr, warn.conflicts = FALSE, quietly = TRUE)
library(dplyr, warn.conflicts = FALSE, quietly = TRUE)
library(ggplot2)
library(ggpubr)
library(sjPlot)
library(kableExtra)
library(xtable)
library(plotly)
library(ggthemes)
library(ggbeeswarm)

# linear model
library(lme4)
library(afex) # mixed effect model
library(permuco)
library(remotes)
library(permutes)
library(coin)

rm(list=ls())

set.seed(0)

Load Data

Load Maximum LogLikelihood data

d <- read.csv("./max_loglikelihood.csv", sep = ",") %>%
  select(HCPID, best_model, maxLL_diff, 
         Punishment_MostlyPunishment.subj, 
         Reward_MostlyPunishment.subj, 
         Reward_MostlyReward.subj, 
         Punishment_MostlyReward.subj) %>%
  pivot_longer(Punishment_MostlyPunishment.subj:Punishment_MostlyReward.subj,
               names_to = "Condition",
               values_to = "Pswitch") %>% 
  mutate(
    block_type = case_when(
      Condition == "Reward_MostlyPunishment.subj" ~ "Loss",
      Condition == "Punishment_MostlyPunishment.subj" ~ "Loss",
      Condition == "Reward_MostlyReward.subj" ~ "Reward",
      Condition == "Punishment_MostlyReward.subj" ~ "Reward",
      TRUE ~ "NA"), 
    trial_type = case_when(
      Condition == "Punishment_MostlyReward.subj" ~ "Loss",
      Condition == "Punishment_MostlyPunishment.subj" ~ "Loss",
      Condition == "Reward_MostlyReward.subj" ~ "Reward",
      Condition == "Reward_MostlyPunishment.subj" ~ "Reward",
      TRUE ~ "NA"), 
    ) 

#d.trial = d %>% 
#  group_by(HCPID, trial_type, best_model) %>%
#  summarise(Pswitch = mean(Pswitch))

Stat Analysis

Mixed GLM in R

We could fit data using glmer

  • Random effects: trial_type, block_type
  • Fixed effects: best_model
  • id: HCPID
  • binomial distribution

Fitting data to mixed effect linear model, best_model label is significant

summary(m <- glmer(data=d, Pswitch ~ d$best_model  + (trial_type+block_type|HCPID), 
      family = "binomial"))

Let’s try other packages.. apex library

Significant fixed effects:

  • best_model

  • trial_type

2 way interaction effects:

  • best_model : block_type

  • best_model : trial_type

3 way interaction effects:

  • best_model:block_type:trial_type
#apex library
mixed_anova <- aov_ez(
  id = "HCPID",
  dv = "Pswitch",
  data = d,
  between = "best_model",
  within = c("block_type", "trial_type")
)
print(mixed_anova)
## Anova Table (Type 3 tests)
## 
## Response: Pswitch
##                             Effect     df  MSE         F   ges p.value
## 1                       best_model 1, 197 0.10 43.55 ***  .099   <.001
## 2                       block_type 1, 197 0.02      1.22 <.001    .270
## 3            best_model:block_type 1, 197 0.02   6.92 **  .004    .009
## 4                       trial_type 1, 197 0.05   7.33 **  .009    .007
## 5            best_model:trial_type 1, 197 0.05 41.10 ***  .046   <.001
## 6            block_type:trial_type 1, 197 0.03      0.43 <.001    .515
## 7 best_model:block_type:trial_type 1, 197 0.03 12.16 ***  .009   <.001
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Visualize mixed effects

ggplotly(afex_plot(mixed_anova, x = "best_model", trace = c("trial_type"), #error = "none",
                   mapping = c("color", "shape"), data_arg = list(color = "gray", shape = 1),
                   line_arg = list(color = c('tomato2', 'tomato2', 'forestgreen', 'forestgreen'), size=2),
                   data_plot = T, point_arg = list(size = 5, 
                                                   color = rep(c('tomato2', 'forestgreen'), each=2), 
                                                   shape = rep(c('x', '.'), each=2))) +
           labs(title = "Behaviror: Trial effect between two models", 
                x = "ACT-R Model Type", y = "Probability of Switch") + 
           theme_bw())
ggplotly(afex_plot(mixed_anova, x = "best_model", trace = c("block_type"), #error = "none",
                   mapping = c("color", "shape"), data_arg = list(color = "gray", shape = 1),
                   line_arg = list(color = c('tomato2', 'tomato2', 'forestgreen', 'forestgreen'), size=2),
                   data_plot = T, point_arg = list(size = 5, 
                                                   color = rep(c('tomato2', 'forestgreen'), each=2), 
                                                   shape = rep(c('x', '.'), each=2))) +
           labs(title = "Behaviror: Block effect between two models", 
                x = "ACT-R Model Type", y = "Probability of Switch") + 
           theme_bw())

Similarly, fit data into mixed effect model

m1 <- mixed(Pswitch ~ best_model + (trial_type + block_type | HCPID), d, 
            method = "LRT", 
            all_fit = TRUE, 
            family = binomial)
nice(m1)
## Mixed Model Anova Table (Type 3 tests, LRT-method)
## 
## Model: Pswitch ~ best_model + (trial_type + block_type | HCPID)
## Data: d
## Df full model: 8
##       Effect df     Chisq p.value
## 1 best_model  1 63.59 ***   <.001
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
xtable(anova(m1))
## % latex table generated in R 4.0.3 by xtable 1.8-4 package
## % Fri May 13 22:04:54 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{lrrrr}
##   \hline
##  & Df & Chisq & Chi Df & Pr($>$Chisq) \\ 
##   \hline
## best\_model & 7 & 63.59 & 1 & 0.0000 \\ 
##    \hline
## \end{tabular}
## \end{table}
summary(m1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Pswitch ~ best_model + (trial_type + block_type | HCPID)
##    Data: data
## Control: ctrl
## 
##      AIC      BIC   logLik deviance df.resid 
##   1027.5   1064.9   -505.7   1011.5      788 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.20028 -0.22055  0.05622  0.42645  1.52324 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev.  Corr       
##  HCPID  (Intercept) 0.000e+00 0.000e+00            
##         trial_type1 5.095e-15 7.138e-08   NaN      
##         block_type1 1.981e-14 1.407e-07   NaN -0.11
## Number of obs: 796, groups:  HCPID, 199
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.23828    0.07699  -3.095  0.00197 ** 
## best_model1 -0.60340    0.07699  -7.837  4.6e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## best_model1 -0.211
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

Permutation Test

Next, we could use apex package to perform permutation test.

d.block = d %>% 
  group_by(HCPID, best_model, block_type) %>%
  summarise(Pswitch = mean(Pswitch), .groups = "drop") %>%
  spread(key = block_type, value = Pswitch) %>%
  mutate(block_effect = Reward - Loss) %>%
  select(HCPID, best_model, block_effect) %>%
  spread(key = best_model, value = block_effect)

d.trial = d %>% 
  group_by(HCPID, best_model, trial_type) %>%
  summarise(Pswitch = mean(Pswitch), .groups = "drop") %>%
  spread(key = trial_type, value = Pswitch) %>%
  mutate(trial_type = Reward - Loss) %>%
  select(HCPID, best_model, trial_type) %>%
  spread(key = best_model, value = trial_type)

d.block.m1 <- d.block %>% select(m1) %>% drop_na()
d.block.m2 <- d.block %>% select(m2) %>% drop_na()

d.trial.m1 <- d.trial %>% select(m1) %>% drop_na()
d.trial.m2 <- d.trial %>% select(m2) %>% drop_na()
compare.2.vectors(d.block.m1$m1, d.block.m2$m2, paired = FALSE, na.rm = FALSE,
  tests = c("parametric", "nonparametric"), coin = TRUE,
  alternative = "two.sided",
  wilcox.exact = NULL, wilcox.correct = TRUE)
## $parametric
##    test test.statistic test.value  test.df           p
## 1     t              t  -2.630521 197.0000 0.009199408
## 2 Welch              t  -2.661803 152.8869 0.008603209
## 
## $nonparametric
##              test test.statistic  test.value test.df          p
## 1 stats::Wilcoxon              W 3604.000000      NA 0.01318405
## 2     permutation              Z   -2.592059      NA 0.00959000
## 3  coin::Wilcoxon              Z   -2.480040      NA 0.01244000
## 4          median              Z   -2.407694      NA 0.01776000
compare.2.vectors(d.trial.m1$m1, d.trial.m2$m2, paired = FALSE, na.rm = FALSE,
  tests = c("parametric", "nonparametric"), coin = TRUE,
  alternative = "two.sided",
  wilcox.exact = NULL, wilcox.correct = TRUE)
## $parametric
##    test test.statistic test.value  test.df            p
## 1     t              t   6.411169 197.0000 1.045368e-09
## 2 Welch              t   6.534571 156.1739 8.553576e-10
## 
## $nonparametric
##              test test.statistic  test.value test.df           p
## 1 stats::Wilcoxon              W 6826.500000      NA 7.72669e-09
## 2     permutation              Z    5.846384      NA 0.00000e+00
## 3  coin::Wilcoxon              Z    5.775596      NA 0.00000e+00
## 4          median              Z    4.338846      NA 2.00000e-05
perm.lmer(Pswitch ~ best_model + (trial_type + block_type | HCPID), data=d, family = binomial(),  type = "regression")
##         Factor      LRT       beta         z  p
## 1  (Intercept) 187.3028 -0.4147827 -4.575609 NA
## 2 best_modelm2 115.7080  0.6174496  4.139432  0

We could also simulate mean difference and calculate the p-value by creating our own permutation test

Block effect: Pswitch(Reward Block) - Pswitch(Loss Block)

Trial effect: Pswitch(Reward Trial) - Pswitch(Loss Trial)

Mean Difference of Block Type: Block effect(M1) - Block effect(M2)

Mean Difference of Trial Type: Trial effect(M1) - Trial effect(M2)

shuffle_labels <- function(d) {
  res <- d %>% 
    select(HCPID, best_model) %>% 
    unique()
  res$best_model <- sample(res$best_model, replace = T)
  res <- d %>% select(-best_model) %>%
    left_join(res, by = "HCPID")
  return (res)
}

block_effect <- function(dat) {
  res = dat %>% group_by(best_model, block_type) %>%
    summarise(Pswitch = mean(Pswitch), .groups = 'drop') %>%
    pivot_wider(names_from = block_type, values_from = Pswitch) %>%
    mutate(block_effect.RL = Reward-Loss) %>%
    select(best_model, block_effect.RL) %>%
    pivot_wider(names_from = best_model, values_from = block_effect.RL) %>%
    mutate(diff.m1m2 = m1-m2) 
  return (res)
}

trial_effect <- function(d) {
  res = d %>% group_by(best_model, trial_type) %>%
    summarise(Pswitch = mean(Pswitch), .groups = 'drop') %>%
    pivot_wider(names_from = trial_type, values_from = Pswitch) %>%
    mutate(trial_effect.RL = Reward-Loss) %>%
    select(best_model, trial_effect.RL) %>%
    pivot_wider(names_from = best_model, values_from = trial_effect.RL) %>%
    mutate(diff.m1m2 = m1-m2) 
  return (res)
}
my_block_effect <- block_effect(d)
my_trial_effect <- trial_effect(d)

# simulate random block effect
n <- 5000
simulated_block_effect = c()
simulated_trial_effect = c()

for (x in 1:n) {
  d0 <- shuffle_labels(d)
  curr.block <- block_effect(d0)
  curr.trial <- trial_effect(d0)
  simulated_block_effect <- c(simulated_block_effect, curr.block$diff.m1m2)
  simulated_trial_effect <- c(simulated_trial_effect, curr.trial$diff.m1m2)
}

Calculate one-side and two-side p-value

# Let's compute the p-value
if (my_block_effect$diff.m1m2 < 0) {
  pvalue.block <- sum(my_block_effect$diff.m1m2 > simulated_block_effect) / n
} else {
  pvalue.block <- sum(my_block_effect$diff.m1m2 <= simulated_block_effect) / n
}
paste("Block Type: The one-sided p-value is ", round(pvalue.block,3))
## [1] "Block Type: The one-sided p-value is  0.005"
paste("Block Type: The two-sided p-value is ", 2 * round(pvalue.block,3))
## [1] "Block Type: The two-sided p-value is  0.01"
if (my_trial_effect$diff.m1m2 < 0) {
  pvalue.trial <- sum(my_trial_effect$diff.m1m2 > simulated_trial_effect) / n
} else {
  pvalue.trial <- sum(my_trial_effect$diff.m1m2 <= simulated_trial_effect) / n
}
paste("Trial Type: The one-sided p-value is ", round(pvalue.trial,3))
## [1] "Trial Type: The one-sided p-value is  0"
paste("Trial Type: The two-sided p-value is ", 2 * round(pvalue.trial,3))
## [1] "Trial Type: The two-sided p-value is  0"
#pvalue.greater <- sum(my_block_effect$diff.m1m2 > simulated_block_effect) / n
#pvalue.smaller <- sum(my_block_effect$diff.m1m2 <= simulated_block_effect) / n

#paste("One-sided p-value (greater) is", pvalue.greater)
#paste("One-sided p-value (smaller) is", pvalue.smaller)

Plot the distribution of mean difference of Block Type between declarative model vs. procedural model

gghistogram(simulated_block_effect, add = "mean", add_density = F, bins = 50, rug = T, fill = 'gray',
            title = "Histogram of simulated mean difference: Block Type", xlab = "mean_difference", ylab="simulated counts") +
  geom_vline(aes(xintercept=my_block_effect$diff.m1m2), linetype="dashed", color = "red") +
  annotate(geom = "text", x = my_block_effect$diff.m1m2*0.8, y = 200, 
           label = paste("observed p = ", 2 * round(pvalue.block,3)))

Plot the distribution of mean difference of Trial Type between declarative model vs. procedural model

gghistogram(simulated_trial_effect, add = "mean", add_density = F, bins = 50, rug = T, fill = 'gray',
            title = "Histogram of simulated mean difference: Trial Type", xlab = "mean_difference", ylab="simulated counts") +
  geom_vline(aes(xintercept=my_trial_effect$diff.m1m2), linetype="dashed", color = "red") +
  annotate(geom = "text", x = my_trial_effect$diff.m1m2*0.8, y = 200, 
           label = paste("observed p = ", 2 * round(pvalue.trial,3)))

Compare different ways of calculating MaxLL

There are different ways of calculating maxLL. Next we could explore which way gives us most distinguishable model labels

  • Group by Block Type and Trial Type

  • Group by Block Type only

  • Group by Trial Type only

d.block <- read.csv("./max_loglikelihood_block.csv", sep = ",") %>%
  select(HCPID, best_model, maxLL_diff, 
         MostlyPunishment.subj, 
         MostlyReward.subj) %>%
  pivot_longer(MostlyPunishment.subj:MostlyPunishment.subj,
               names_to = "Condition",
               values_to = "Pswitch") %>% 
  mutate(
    block_type = case_when(
      Condition == "MostlyReward.subj" ~ "Reward",
      Condition == "MostlyPunishment.subj" ~ "Loss",
      TRUE ~ "NA")) 

d.trial <- read.csv("./max_loglikelihood_trial.csv", sep = ",") %>%
  select(HCPID, best_model, maxLL_diff, 
         Reward.subj, 
         Punishment.subj) %>%
  pivot_longer(Reward.subj:Punishment.subj,
               names_to = "Condition",
               values_to = "Pswitch") %>% 
  mutate(
    block_type = case_when(
      Condition == "Reward.subj" ~ "Reward",
      Condition == "Punishment.subj" ~ "Loss",
      TRUE ~ "NA"), 
    trial_type = case_when(
      Condition == "Reward.subj" ~ "Reward",
      Condition == "Punishment.subj" ~ "Loss",
      TRUE ~ "NA")) 

Below plot shows the differences of MaxLL between two models. Different color indicates different ways of calculating MaxLL.

  • For Both, the maxLL was calcualted with aggregated Pswitch by both block type and trrial type

  • For block_type, the maxLL was calcualted with aggregated Pswitch only by block type

  • For trial_type, the maxLL was calcualted with aggregated Pswitch only by trial type

x-axis is the absolute difference between maxLL of model1 (declarative model) and maxLL of model2 (procedural model).

Take-away: If only look at block effect (blue density curve), two model are more distinguisable from each other.

maxLL_diff = data.frame(maxLL_type=c(rep("both", length(d$maxLL_diff)), 
                               rep("block_type", length(d.block$maxLL_diff)), 
                               rep("trial_type", length(d.trial$maxLL_diff))), 
                        maxLL_diff = c(d$maxLL_diff, 
                                       d.block$maxLL_diff, 
                                       d.trial$maxLL_diff)) %>%
  mutate(value_type=factor(maxLL_type, levels = c("both", "block_type", "trial_type")))
 
ggdensity(data=maxLL_diff, x="maxLL_diff", fill = "maxLL_type", palette = "Set1", add = "mean", 
          color = "darkgray", alpha = .5,
          xlab = "max loglikelihood abs(model1 - model2)",
          title = "Density Plot of Max LogLikelihood Differences") +
  theme_pander()

Let’s check whether three different ways give us reliable individual strategy selection.

best_models = d.trial %>% select(HCPID, best_model) %>% unique() %>% 
  left_join(d.block %>% select(HCPID, best_model) %>% unique(), by = "HCPID", suffix = c(".trial", ".block")) %>%
  left_join(d %>% select(HCPID, best_model) %>% unique(), by = "HCPID") %>%
  mutate(m1_count = if_else(best_model=='m1', 1, 0) + 
           if_else(best_model.trial=='m1', 1, 0) + 
           if_else(best_model.block=='m1', 1, 0)) %>%
  arrange(m1_count)

Unfortunately, if we use block type or trial type solely to calcualte maxLL, it simply believes all subjects use m1. This is clearly not what we wants.

best_models %>%
  count(best_model)
## # A tibble: 2 x 2
##   best_model     n
## * <chr>      <int>
## 1 m1           126
## 2 m2            72
best_models %>%
  count(best_model.trial)
## # A tibble: 2 x 2
##   best_model.trial     n
## * <chr>            <int>
## 1 m1                 196
## 2 m2                   2
best_models %>%
  count(best_model.block)
## # A tibble: 1 x 2
##   best_model.block     n
## * <chr>            <int>
## 1 m1                 198

Apply weights to labels

d %>% select(HCPID, best_model, maxLL_diff) %>%
  arrange(best_model, desc(maxLL_diff))
## # A tibble: 796 x 3
##    HCPID       best_model maxLL_diff
##    <chr>       <chr>           <dbl>
##  1 164131_fnca m1               7.50
##  2 164131_fnca m1               7.50
##  3 164131_fnca m1               7.50
##  4 164131_fnca m1               7.50
##  5 169141_fnca m1               5.53
##  6 169141_fnca m1               5.53
##  7 169141_fnca m1               5.53
##  8 169141_fnca m1               5.53
##  9 155231_fnca m1               5.16
## 10 155231_fnca m1               5.16
## # … with 786 more rows